function [acc] = heateqcn(x0, xn, dx, t0, tn, dt)

% [acc] = heateqfinitediff(x0, xn, dx, t0, tn, dt)
%
%   Solves the heat equation diffusion problem in 1d using finite
%   differences on the domain x = [x0, xn] from t = [t0, tn].
%
%               (-4*(x-5)^2)
%   F(x, 0) = e^
%
%   fx0 = initial temperature function
%   dx = spatial step
%   dt = time step

% 0. Setup.
syms fx0x;
fx0 = inline('exp(-4*(x-5)^2)');

nj = 1 + (tn - t0)/dt; % (j) y-axis
ni = 1 + (xn - x0)/dx; % (i) x-axis

acc = zeros(nj+1, ni+1);

% 1. Calculate 'r' for coeff. matrix.
r = dt/(dx^2);

% 2. Create coeff. matrix.
n = ni - 2;
A = zeros(n, n);

for k = 1:n
    if k < n
        A(k, k+1) = -r/2;
        A(k+1, k) = -r/2;
    end
    
    A(k, k) = 1+r;
end

% 3. Create b vector.
b = zeros(n, 1);

% 4. Setup solution matrix.
for i = 1:ni
    acc(1, i) = feval(fx0, dx*(i-1));
    acc(nj+1, i) = dx*(i-1);
end

for j = 1:nj
    acc(j, 1) = 0;
    acc(j, ni) = 0;
    acc(j, ni+1) = dt*(j-1);
end

% 4. Begin iterating.
for j = 2:nj
    
    % Copy previous timestep to current b vector.
    b = zeros(n, 1);
    for k = 1:n
        i = k+1;
        
        if k > 1
            b(k) = b(k) + r/2*acc(j-1, i-1);
        end
        
        b(k) = b(k) + (1-r)*acc(j-1, i);
        
        if k < n
            b(k) = b(k) + r/2*acc(j-1, i+1);
        end
    end
    
    acc(j, 2:ni-1) = linsolve(A, b);
end